STOCHASTIC AND DETERMINISTIC MOLECULAR DYNAMICS DERIVED FROM 
THE TIME-INDEPENDENT SCHRODINGER EQUATION 



ANDERS SZEPESSY 

Abstract. Ehrenfest, Born-Oppenheimer, Langcvin and Smoluchowski dynamics are shown to be accurate 
approximations of time-independent Schrodinger observables for a molecular system avoiding caustics, in 
the limit of large ratio of nuclei and electron masses, without assuming that the nuclei are localized to 
vanishing domains. The derivation, based on a Hamiltonian system interpretation of the Schrodinger equation 
and stability of the corresponding Hamilton-Jacobi equation, bypasses the usual separation of nuclei and 
electron wave functions, includes crossing electron eigenvalues, and gives a different perspective on the Born- 
Oppenheimer approximation, Schrodinger Hamiltonian systems, stochastic electron equilibrium states and 
numerical simulation in molecular dynamics modeling. 
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1. The Schrodinger and molecular dynamics models 
The time-independent Schrodinger equation 

(1.1) H{x,X)<S>{x,X) ^ E<i>{x,X), 

models nuclei-electron systems and is obtained from minimization of the energy in the solution space of 
wave functions, cf. [351 S?! |31 [301 [TO]- It is an eigenvalue problem for the energy E E H oi the system in 
the solution space, described by wave functions, $ : R"^"^ x R'^^ — ^ C, depending on electron coordinates 
X = {x^, . . . , x"^) £ R'^'-', nuclei coordinates X = {X^ , . . . , X^) G R'^^, and a Hamiltonian operator H{x, X) 

1 ^ 

(1.2) H{x,X) = Vix,X)--M-^Y.^X"- 

n=l 

The nuclei masses M are assumed to be large and the interaction potential V, independent of Af , is in the 
canonical setting (neglecting relativistic and magnetic effect), 

f. o\ y{x,X)^ ~^Y.j=l^xi +Y.l<k<j<J Ixk-xJ] 

composed of the kinetic energy of the electrons, the electron-electron repulsion, the electron-nuclei attraction, 
and the repulsion of nuclei (with charge ^n), in the Hartree atomic units where the electron mass, electron 
charge, reduced Planck constant, and the Coulomb force constant (47reo)~^ all are one. The mass of the 
nuclei, which are much greater than one (electron mass), are the diagonal elements in the diagonal matrix 
M. 

An essential feature of the partial differential equation ()l.ip is the high computational complexity to 
determine the solution, in an antisymmetric/symmetric subset of the Sobolev space H^{'R^^^^^^^'>). Wave 
functions depend also on discrete spin states 

(1.4) $(a;\CTi,...,x"',CTj,Xi,I]i,...,XAr,SAr), 

which effect the solutions space: each electron spin aj can take two different values and each nucleus can be 
in a finite set of spin states I]„; the Pauli exclusion principle restricts the solutions space to wave functions 
satisfying the antisymmetry/symmetry 

$(. . . ,x-',crj, . . . ,a;'',crfc, . . .) = -$(... ,2;'', CTfc, ... ,a;^crj, .. .) for any 1 < j, fc < J 

and for any pair of nuclei n and m, with A nucleons and the same number of protons and neutrons, 

$(..., X"\ E™, . . . , X", an, . . .) - (-1)^<I>(. . . , X", E„, . . . , X'", S™, . . .) 

cf. [To]. We simplify the notation by writing $(a;, X) instead of the more complete ()1.4|) . since the Hamiltonian 
H does not depend on the spin. The time-independent Schrodinger equation has convincing agreement with 
experimental results, as the basis for computational chemistry and solid state physics. An attractive property 
of the Schrodinger equation (jl.ip is the precise definition of the Hamiltonian and the solutions space, without 
unknown parameters. The agreement with measurements can be further improved by including relativistic 
and magnetic effects, cf. [TO] . 

In contrast to the Schrodinger equation, a molecular dynamics model of nuclei X : [0,r] — R'^^, with 
a given potential Vp : R'^^ — > R, can be computationally studied also for large N by solving the ordinary 
differential equation 

(1.5) MXr = -dxVp{Xr). 

This computational and conceptual simplification motivates the study to determine the potential and its 
implied accuracy by a derivation of molecular dynamics from the Schrodinger equation, as started already in 
the 1920's with the seminal Born-Oppenheimer approximation [6]. The purpose here is to contribute to the 
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current understanding of such derivations, by showing convergence rates under new assumptions. The precise 
aim in this paper is to estimate the error 

(1.6) [ g{X)<i>{x, Xy^x, X)dxdX - \im T^^ [ g{Xr)dT 

of a given position observable / g{X)(^{x^X)*(^{x,X)dxdX of the time-indepedent Schrodinger equation 

(jl.l[) approximated by the corresponding molecular dynamics observable Wuyt^oo g{XT)dT, which is 

computationally cheaper to evaluate for several nuclei. 

A useful sub step to derive molecular dynamics from the Schrodinger equation is Ehrenfest dynamics, for 
classical ah initio motion of the nuclei coupled to Schrodinger dynamics for the electrons, 

(-^ 7) MXl' = - /j^3,7 rAx, Xr) dx^V{x, Xr) ^r) dx 

with the initial normalization J-^^aj 4>q{^^ ^o)4>o{x, XQ)dx — 1. The Ehrenfest dynamics (II. 7p has been derived 
from the time-dependent Schrodinger equation through the self consistent field equations, see [71 |40l [51]. 
Equation ()1.7p can be used for ab initio computation of molecular dynamics, cf. [40( 137] . A next step is 
the zero-order Born-Oppenheimer approximation, where Xt solves the classical ah initio molecular dynamics 
(|1.5p with the potential Vp : R^^ — ^ R determined as an eigenvalue of the electron Hamiltonian V{-,X) for 
a given nuclei position X, that is = Aq and V{-,X)il;o{-,X) = \o{X)i/jo{-,X) for an electron eigenfunction 
ipo{',X) e i^(R'^'^), for instance the ground state. 

The Born-Oppenheimer expansion [6] is an approximation of the solution to the time-independent 
Schrodinger equation which is shown |27[ I32j to solve the time-independent Schrodinger equation approxi- 
mately. This expansion, analyzed by the methods of multiple scales, pseudo differential operators and spectral 
analysis in [271 1321 US] , can be used to study the approximation error (|1.6|) , cf . Section 12.4.11 Although in 
the literature it seems easier to find precise statements on the error in observables for the other setting of 
the time-dependent Schrodinger equation. Instead of an asymptotic expansion, we use a different method 
based on a Hamiltonian dynamics formulation of the time-independent Schrodinger eigenfunction and the 
stability of the corresponding perturbed Hamilton- Jacobi equations viewed as a hitting problem. Another 
motivation for our method is that it forms a sub step in trying to estimate the approximation error using 
only information available in molecular dynamics simulations. 

The related problem of approximating observables to the time-dependent Schrodinger by the Born Oppen- 
heimer expansions is well studied, theoretically [8ll43j and computationally [35^ using the Egorov theorem. The 
Egorov theorem shows that finite time observables of the time-dependent Schrodinger equation are approx- 
imated with O(Af-i) accuracy by the zero-order Born-Oppenheimer dynamics, with an electron eigenvalue 
gap. In the special case of a position observable and no electrons (i.e. V=V(X) in (|1.2p ). the Egorov theorem 
states that 



/ giXMX,tr^X,t)dX- f 



(1.8) / g{XMX,ty^X,t)dX- / g(Xt)$(Xo,0)*$(Xo,0)dXo 



< CfAr 



where $(X, t) is a solution to the time-dependent Schrodinger equation iM ^/^9t$(-,t) — H^{-,t) with 
the Hamiltonian (|1.2p and the path Xt is the nuclei coordinates for the dynamics with the Hamiltonian 
|Xp/2 + V{X). If the initial wave function ^{X, 0) is the eigenfunction in (jl.ip the first term in (|1.8p reduces 
to the first term in (II. 6|) and the second terms can also become the same in an ergodic limit; however since 
we do not know that the parameter Cf (bounding an integral over {0,t)) is bounded for all time we cannot 
directly conclude an estimate for (|1.6p from ()1.8|) . In this perspective, to study the time- independent instead 
of the time-dependent Schrodinger equation has the important differences that 

• the infinite time study of the Born-Oppenheimer dynamics can be reduced to a finite time hitting 
problem. 
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• the computational and theoretical problem of specifying initial data for the Schrodinger equation is 
avoided, and 

• computational cheap evaluation of a position observable g{X) is possible using the time average 
limT^oo Jq 9i^T)dT/T along the solution path Xr- 

This paper derives the Ehrenfest dynamics (jl.7|) and the Born-Oppenheimer approximation from the time- 
independent Schrodinger equation and establishes convergence rates for molecular dynamics approxi- 
mations to time- independent Schrodinger observables under simple assumptions excluding so called caustic 
points, where the Jacobian determinant det dXt/dXo of the Eulerian-Lagrangian transformation of X— paths 
vanish. As mentioned, the main new analysis idea is to write the time-independent Schrodinger equation 
as a Hamiltonian system and analyze the approximations by comparing their Hamiltonians with the 
Schrodinger Hamiltonian, using the theory of Hamilton- Jacobi partial differential equations; the problematic 
infinite time evolution of perturbations in the dynamics is solved by viewing it as a finite time hitting problem 
for the Hamilton- Jacobi equation. Another difference is we analyze the transport equation as a time-dependent 
Schrodinger equation in contrast to the traditional rigorous and formal asymptotic expansions. 

The main inspiration for this paper is [131 [H [S] and the semi-classical WKB analysis in 01] • the work 
|42[ m [5] derives the time-dependent Schrodinger dynamics of an x— system, = Hi^!, from the time- 
independent Schrodinger equation (with the Hamiltonian Hi{x) + 6H{x, X)) by a classical limit for the 
environment variable A", as the coupling parameter 5 vanishes and the mass M tends to infinity; in particular 
[m m [S] show that the time derivative enters through the coupling of with the classical velocity. Here we 
refine the use of characteristics to study classical ab initio molecular dynamics, where the coupling does not 
vanish, and we establish error estimates for Born-Oppenheimer, Ehrenfest, surface-hopping and stochastic 
approximations of the Schrodinger observables in the case of no caustics present. The small scale, introduced 
by the perturbation 

-(2M)-i^Ax„ 

n 

of the potential V , is identified in a modified WKB eikonal equation and analyzed trough the corresponding 
transport equation as a time-dependent Schrodinger equation along the eikonal characteristics. This modified 
WKB formulation reduces to the standard semi-classical approximation, cf. [41 , for the case of a potential 
function V{X) g R (depending only on nuclei coordinates) but becomes different in the case of operator 
valued potentials studied here. The global analysis of WKB functions started by Maslov in the I960' [41] 
and lead to the subject Geometry of Quantization, relating global classical paths to eigenfunctions of the 
Schrodinger equation, cf. jl4j . The analysis here, based on a Hamiltonian system interpretation of the time- 
independent Schrodinger equation and stability of the corresponding Hamilton- Jacobi equation, bypasses the 
usual separation of nuclei and electron wave functions in the time-dependent self consistent field equations 

[niioKsi]. 

Theorems 16.11 and 16.21 show that observables based on a single WKB Schrodinger eigenstate are approx- 
imated by observables from the Ehrenfest dynamics and the zero-order Born-Oppenheimer dynamics with 
error 0{M~"), using that these approximate solutions generate approximate eigenstates to the Schrodinger 
equation: in the case of no caustics and the electron eigenvalues satisfying a spectral gap condition, there 
holds a = 1 — 6 for any S > 0; for non degenerate electron eigenvalue crossings the convergence rate is reduced 
to a = 1/2 — (5 for Born-Oppenheimer dynamics and a = 3/4 for Ehrenfest dynamics, using the stationary 
phase method in Section 17.41 The results are based on the Hamiltonian (II. 2p with any potential V that is 
smooth in A, e.g. a regularized version of the Coulomb potential (11.31) . The derivation does not assume 
that the nuclei are supported on small domains; in contrast, derivations based on the time-dependent self 
consistent field equations require nuclei to be supported on small domains. The reason that small support is 
not needed here comes from the combination of the characteristics and sampling from an equilibrium density, 
that is, the nuclei paths behave classically although they may not be supported on small domains, in the case 
of no caustics. Section [5] shows that caustics couple the WKB modes, as is well known from geometric optics. 
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see [3T1|4T], and generate non orthogonal WKB modes that are coupled in the Schrodmger density; otherwise, 
without caustics the Schrodinger density is asymptotically decoupled into a simple sum of individual WKB 
densities. Remark 16.31 relates the approximation results to the accuracy of symplectic numerical methods for 
molecular dynamics. 

Section [3] shows that the Ehrenfest dynamics is formally the same when derived from the time-independent 
and time-dependent Schrodinger equations; a unique property of the time-independent Schrodinger equation 
we use is the interpretation of its dynamics Xt G returning to a co-dimension one surface / and thereby 
reducing the dynamics to a hitting time problem with finite time excursions from /. Another advantage with 
molecular dynamics approximating an eigenvalue, is that stochastic perturbations of the electron ground 
state can be interpreted as a Gibbs distribution of degenerate nuclei-electron eigenstates of the Schrodinger 
eigenvalue problem The time- independent eigenvalue setting also avoids the issue on "wave function 

collapse" to an eigenstate, present in the time-dependent Schrodinger equation. 

The model (|1.5|) simulates dynamics at constant energy M\X\'^/2 + Vp{X), constant number of particles 
N and constant volume, i.e. the microcanonical ensemble. The alternative to simulate with constant number 
of particles, constant volume and constant temperature T, i.e. the canonical ensemble, is possible for instance 
with the stochastic Langevin dynamics 

. , dXr = Vrdr 

^ ' Mdvr ^-dxVp{Xt)dT-KvtdT+{2TKf/^dWr, 

where Wr is the standard Brownian process (at time r) in with independent components and X is a 
positive friction parameter. When the observable only depends on the nuclei positions, i.e. not on the nuclei 
velocities or the correlation of positions at different times, the Smoluchowski dynamics 

(1.10) dXr = -dxVpiXr) + {2Tfl'^dWr 

is a simplified alternative to Langevin dynamics, cf. [9]. 

Section 17.21 shows that the Schrodinger eigenvalue problem can be written as a Hamiltonian molecular 
dynamics system and that the equilibrium Gibbs distribution of the Ehrenfest Hamiltonian system approxi- 
mates the Gibbs distribution of the Hamiltonian Schrodinger dynamics with accuracy 0{M~^) for a spectral 
gap case. Theorem 18.31 shows that Langevin and Smoluchowski dynamics accurately approximate observ- 
ables of Schrodinger and Ehrenfest Gibbs equilibrium dynamics, when the observable depends only on the 
nuclei positions but not their correlation at different time. The derivation uses a classical equilibrium Gibbs 
distribution, based on a probability distribution to be in any WKB eigenstate (corresponding to a degener- 
ate eigenvalue of the Schrodinger equation (jl.ip ). randomly perturbed from the electron ground state; the 
derivation uses an other assumption of a spectral gap and no caustics. The main idea in the theorem is the 
formulation of a classical Gibbs equilibrium distribution of eigenstates, motivated by nuclei acting as heat 
bath for the electrons in the quantum Ehrenfest and Schrodinger Hamiltonian systems. Theorem 18.51 shows 
that stochastic perturbations of the ground state can also generate a temperature dependent contribution to 
the drift, depending on a spectral gap of the electron eigenvalues. 

I hope that these ideas can be further developed to better understand molecular dynamics simulations, for 
instance 

- I assume that the electron eigenfunction is smooth enough as a function of X, which holds when 
the potential V is smooth in X - to find a possibly reduced convergence rate in the Coulomb case 
()1.3|) will require a more careful study, 

- it would be desirable to understand the effect of caustics, as it is for systems without electrons (or 
with electrons as heavy as the nuclei) in the so called semi-classical limit, cf. [H]. 

We use the notation ■0(x,X) — 0{M~°-) also for complex valued functions, meaning that |?/'(x,X)| — 
0{M^") holds uniformly in x and X. 
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2. Ehrenfest dynamics derived from the time-independent Schrodinger equation 

2.1. Exact Schrodinger dynamics. Assume for simplicity that all nuclei have the same massQ The sin- 
gular perturbation — (2A/)^^ of the potential V introduces an additional small scale A/^^/^ of high 
frequency oscillations, as shown by a WKB-expansion, see [45) [29l [28l [13] . We will construct solutions to 
IlT]) in such WKB-form 

(2.1) ^x,X)^ijix,X)e''""'"'>^''\ 

where the wave function is complex valued, the phase 9 is real valued and the factor Af^/^ is introduced 
to have well defined limits of ip and as Af — oo. The standard WKB-construction [51] [55] is based on 
a series expansion in powers of A/^/^ which solves the Schrodinger equation with arbitrary high accuracy. 
We introduce instead of an asymptotic solution an actual solution based on a time-dependent Schrodinger 
transport equation. This transport equation reduces to the formulation in [IT] for the case of a potential 
function V{X) € R (depending only on nuclei coordinates X G R"^^) and modifies it for the case of a self- 
adjoint potential operator V{-,X) on the electron space L^(R'^'-') focused on here; the second difference is 
that we analyze the transport equation as a time-dependent Schrodinger equation instead of as an asymptotic 
expansion. In the next section we use a linear combination of such eigensolutions. The WKB-solution satisfies 
the Schrodinger equation ()l.ip provided that 

= (i/ - S)V^e'*^'''«W 

(2.2) ^[O-^ + V^E]^ 

-237 Ax.V' - ig^T^ j:jidx^i^dx.9+ ^i^dx.x^9))e^^'''''>(^^ . 

We will see that only eigensolutions $ that correspond to dynamics without caustics correspond to such a 
single WKB-mode, as for instance when the eigenvalue E is inside an electron eigenvalue gap. Introduce the 
complex-valued scalar product 



v{x, ■)*w{x, ■)dx = (vlw) 

on L^(R^'') and the notation X»Y for the standard scalar product on R^^ . To find an equation for 9, multiply 
(|2.2p by ■!/'*e~**^ ^ and integrate over R'^'^; similarly take the complex conjugate of (|2.2[) . multiply by 
^giM e(x) integrate over R'^'^; and finally add the two expressions to obtain 

=20-^^^-E) ^-iP + ^P-V^p + Vip-^ 
1 I ^ ^ , 

-217 • (E, Ax. V) + (E, Ax. V-) • V-) 
(2.3) -Jim[ f- jdx^ • 9x9) - {dxi^ . dx9) ■ jj^ 



( ^(^ . V; - ^ . V;) dx. X. e) . 



2Mi 

=0 

The purpose of the phase function 9 is to generate an accurate approximation in the limit as A^ — >■ oo: 
therefore we define 9 by the formal limit of (|2.3p as Af — )■ c», which is the Hamilton- J acobi equation (also 
called the eikonal equation) 

(2.4) ^ =E~V,, 



^If this is not the case, change to new coordinates M^^'^X^ = M^^'^X'', which transform the Hamiltonian to the form we 
,1/2 



want V{x, mI'^M-^/^X) - (2Mi)-i J2n=i ^X' 
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where the function Vq : R'^^ — >■ R is 

This Eikonal equation (|2.4p reduces to the standard formulation 
(2.5) (^^_i5 + y)^ = 0, 

cf. [3T], when the potential function V{X) E R depends only on nuclei coordinates X e R^^^ but (|2.4p is 
different for the case with an operator valued self-adjoint potential V{-,X) on the electron space L^(R'^-^); in 
particular the standard condition (I2.5|) requires ip to be an eigenfunction of V. The reason we use the Eikonal 
equation (|2.4I) instead of (12.51) is that the corresponding transport equation to leading order becomes a variant 
of the Ehrenfest Hamiltonian dynamics (|1.7p . as shown in the remainder of this section. This section is also 
a step towards the more basic construction in Section [721 were we modify Vq in (|7.7|) with an asymptotically 
negligible term, to write the time-independent Schodinger equation as a Hamiltonian system, so that the 
combination of the Eikonal equation and the corresponding transport equation forms a Hamiltonian system. 

For the energy E chosen larger than the potential energy, that is such that E > Vq, the method of 
characteristics, cf. [18) , 



dXt 
dt 



^dx9{Xt) =:pu 



(2.6) ^ =-dxV,{X,), 

^ =\p,\^ = 2{E-Vo{X,)), 

yields a solution (X, p, z) : [0,r] — ^ J7 x R^^ x R to the eikonal equation (|2.4p locally in a neighborhood 
U C R-^^ , for regular compatible data {Xq ,po, zq) given on a 3 — 1 dimensional " inflow" -domain I C U; here 
Zt '■— 9{Xt). Typically the domain / and the data 0\i are not given, unless it is really an inflow domain and 
characteristic paths do not return to J as in a scattering problem. If paths leaving from / return to /, there 
is an additional compatibility of data on /: we have zq = — J* \p^\'^ds + zt, where Xq e / and Xt G /, so that 
Zt = 0{Xt) and pt = dxO{Xt) are determined from zq = 0{Xo) and po — dx6{Xo), and continuing the path to 
subsequent hitting points Xt^ & I, j = 1,2,... determines (^9{Xt-),dxO{Xt-)) from (^9{Xo),dxd{Xo)) . Our 
derivation of approximation error will use such a WKB Ansatz for the approximate Ehrenfest solution, the 
Born-Oppenheimer solution and for the exact Schrodinger solution $. As we shall see, the Ehrenfest and 
Born-Oppenheimer approximations have related (simpler) equations for its characteristics. 

The phase function 9 : U ^ H becomes globally defined in J7 C R^^ when there is a unique characteristic 
path Xt going through each point in U. A globally defined wave function $ can be constructed from a linear 
combination of WKB functions also when caustics are present, using the manifold of phase-space solutions 
{X, p) and Fourier integral operators to relate X and p dependence, see [41] and [14] . We assume in this work 
that characteristics do not collide to avoid multi valued solutions. A simple case without caustics is when the 
potential is such that minxeR(-E' ^ K)(A^)) > 0, in one dimension. Section [4. II presents approximations with 
a linear combination of WKB functions related to so called surface-hopping. 

Definition (j2.4p and equation (|2.2p imply that solves the so called transport equation 

(2.7) _ _i_ ^ A^^^ + _ Vo)iJ - ^ {9x.i^dx.e+ \dx,xoei>). 

3 3 

Time enters into the Schrodinger equation through the characteristics and the chain rule 

OxV • OxO = dxV • —J— = Ti ' 

at dt 

cf. [m m [S] . The second term in the right hand side of (|2.7p can be simplified by the scalar integrating factor 

(2.8) Gt:=e<^S 
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defined along the characteristics from 

(2.9) ^-^E^^^^^^(^*)- 
The integrating factor Gt, defined by (|2.8p and \2.9\ . gives 



Gt dt ■ 

We can also define a function G : Yi?-^ — > R by identifying G{Xt) '■— Gt- This step, with the integrating 
factor G, differs from [42l|4l[5], which approximates the last term in (|2.7p . 

by zero in their case of vanishing coupling between the quantum system and the environment; here the coupling 
between the nuclei and electrons does not vanish. The right hand side in (|2.7p becomes the time derivative 
%M~^I'^G~^ d{G^)/dt and we have derived the time- dependent Schrddinger equation, for the variable tp := 

= {H - E)^ 

= (( - Mh^i' + iV- Vo)^ - ^ GAx. (^G-i))G-i 

^ ^ ^ 

=0 

The density can be written 

(2.11) P — / \M^dx= [ \ij\^dx/G^ 
and therefore the second equation in (|2.6p yields the nuclei dynamics 

if) ■ ip 

The weight function G^ equals the determinant of the first variation dXt / OXq modulo a constant 

(2.12) G^JGl = detidXt/dXo), 

which follows from Liouville's formula, see [H], and in one dimension G^/Gq —pt/po- 

In conclusion, we have rewritten a WKB solution of the time-independent Schrodinger equation in the 
form of the exact Schrddinger dynamics along the characteristics of the Eikonal equation and the transport 
equation 



(2.13) 



X = -dx^ 



The eigenvalue is a parameter in the Hamiltonian for the characteristics of the Hamilton- Jacobi equation 
(|2.4p . In the case when no electrons are present, we have F = Vb, and equation (|2.13p was derived from 
\1.2\ and (|2.ip in [41]. The integrating factor G and its derivative dxG can be determined from (p, dxp, dxxp) 
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along the characteristics by the following characteristic equations obtained from ^-differentiation of p.4p 



dt 



(2.14) 



= -T.j dx-P^dxkp^ - dxrxi'Vo, 



Tfdx^x^p'' = 



Y.jP'dx^x-x-iP^ = T.jP'^X-X'-XiP^ 
= - J2j clx-p-'dxkxiP-' - dx-xip-'dxkp-' - Sx-jc'^x-j V'o, 

and similarly dxxG can be determined from (p, dxP, dxxP, dxxxp)- 

2.2. Approximate Ehrenfest dynamics and densities. We define the approximating Ehrenfest dynamics 
by in ([22]) neglecting the term {2My^ A^j-iA: 

(2.15) ^ + V - E^ti: - j^j:^{dx:i[^dxJ + ^i^dx.x.0), 

and seek, as in p.4p . the approximate phase 6 as the solution to eikonal equation 

\dxO\' 



(2.16) 
where 

Introduce its characteristics 
to rewrite (|2J5)) . as in (fO))) . 



Vn 



dXt 
dt 

dpi 
dt 

dzi 
dt 



^dxkXt) =:pt, 

= -dxVoiXt), 

= \pt? ^2{E-VoiXt)), 



=-{j^i^-{V-VoH)G-' + { 



\dxO\ 



+ Vo~E) ^, 



(2.17) 



for -0 := Gip approximating ip and 



(2.18) 



Gt 



as in (|2.9p (where C is a positive constant for each characteristic). Using that ip ■ ip is conserved (i.e. time- 
independent) in the Ehrenfest dynamics, we can normalize io ip ■ ^ — 1. Note that in the exact dynamics, 
the function ip ■ ^ is not conserved, due to the L^(R^'^) non symmetric source term GAj^j (V'/G) in 

(mSl). We have 



(2.19) 
and 



tp ■ Vip ■0 • Vxp 
ip Pp . ^ 



Vo = — ^ = -0 • ^V- 
tp ■ tp 
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2.3. Comparison of two alternative Ehrenfest formulations. The two different Ehrenfest dynamics 
(l2Tfl) respectively (fTTl) differ in: 

(1) the different time scales t (slow) respectively M^^^t =: r (fast); 

(2) the potentials V — Vq and V in the equations for if; and 0, respectively; and 

(3) the forces dxi'ip ' ^"0) respectively (j) ■ dxV(j) in the momentum equation. 

If -0 solves (|2.17p , the change of variables 

and the property -0 • A-tp — (j) ■ A(j), which holds for observables A not including the nuclei momentum idx (in 
particular for A — V), imply that (j) solves 

(2.20) 

There has been a discussion in the literature [71 [5T] whether the forces should be computed as above in 
((2:201 or as in ([LTl) by 

(221) ^-J^,jr{;Xr)dx^V{;Xr)^i;Xr)dx 

Theorems 16.11 and 16 . 2 1 below show that both formulations p. 201) and (|2.211[T77)) yield accurate approximations 
of the Schrodinger observables, although they are not the same and (I2.21[ II. 7p is closer to the Schrodinger 
observables: the reason that both approximations are accurate is that (|2.2H[T77l) forms a Hamiltonian system 
(as explained below and in Section |6]) which is close the Hamiltonian system for the Schrodinger equation, 
see Theorem 16.11 while the formulation (j2.20l) can be viewed as a Hamiltonian system approximating the 
Born-Oppenheimer dynamics (jl.5l) . using that the wave function ip is close to an electron eigenfunction (see 
Section [7.4p . and the Born-Oppenheimer dynamics approximates Schrodinger observables, see Theorem 16.21 
The formulation (|1.7p and (|2.2ip has the advantage to be a closed Hamiltonian system: the variable 
{X,ipr;p,ipi), with the definition 

d^^O =: (fii, 

and the Ansatz 9 — 9 imply that the Hamilton- Jacobi equation (|2.16p becomes 

^8x9 . 8x9 + (j)r ■ V{X)(^r + 4>^ ' V{X)4>, 

\dx9 • 8x9 + 2-iAf i/Vr • V{X)ipr + 2~^M^/^^i ■ V{X)(fi 
\dx~9 • 8x~9 + ■ V{X)<fr + 2~^M^/^8^J ■ V{X)d^J 

E, 

where the derivative 8^^9{X, (pr) ~ ^i, of the functional 9 : R'^^ x L^{dx) — > R, is the Gateaux derivative in 
L^{dx). Its characteristics form the Hamiltonian system 

Xt ^Pt 

Pt =-^ip{t)-dxV{XtMt) 
Mt) = My^v{Xt)ip,it) 



(2.22) 
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which is the same as the Hamiltonian system (|2.2ip 

Xt ^Pt 

(2.23) pt = ~<pt ■ dxV{Xt)^t 

j^<Pt =V{Xt)cl)t. 

The Hamihonian system yields the equation for the phase 

d = dxO>X + d^J ■iPr=p•p+2(t>^■V(t)^= 2{E - <j>r ■ V(j)r), 

since 9 = 9{X, ipr) is a function of both X and tpr in this formulation. The important property of this 
Hamiltonian dynamics is that -0), with 

solves both the Hamilton- Jacobi equation (j2.16p and ()2.15p . which leads to an approximate solution of the 
Schrodinger eigenvalue problem in (|6.ip - (|6.3p . 

The alternative (|2.20p does not form a closed system, in the sense that the required function dx'ipiXt) is 
not explicitly determined along the characteristics, but can be obtained from values of 0, in a neighborhood 
of Xf, by differentiation. 

2.4. Equations for the density. We note that 

1/2 



/ . \l/2 . 

^ =te) ^ 



shows that the densities p = ip ■ tp (defined in p.lip ) and p := ip ■ ip, in addition to {X,p,tp) and {X,p,tp), 
are needed to determine the wave functions ^ and ip. 

Equation (|2.2I) and the projections in (|2.3p subtracted (instead of added) imply that the density p satisfies 

= J2j dxoipdx^d), 

where denotes the imaginary part of w, and consequently, the density can be determined along a charac- 
teristic using (HHl) and 

p{Xt) =Y.,dx,p{Xt)X^ 

(2.24) - -p{Xt) dx^x.O- M-V2 3(0 . Ax^^) 

= -p{Xt) div-Af-V2^^.cj(^. A;t.0;)p 
- -p(XO^ log - Af- 1/2 cj(^ . Ax.^P). 

Similarly, the Ehrenfest density satisfies the conservation of mass 

O^T.^xApdxJ) 

j 

so that 

(2 25) Z^^^*) ^Ejdx,p{Xt)jc^ 

= -/5(X,)|logG? 

with the solution 

(2.26) p{Xt) = ^, 

G* 
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where C is a positive constant for each characteristic. Note that the derivation of this classical density 
does not need a corresponding WKB equation but uses only the conservation of mass that holds for classical 
paths. The classical density corresponds precisely to the Eulerian-Lagrangian change of coordinates G^/Gg = 
det{dXt/dXo) in (|2.12p . Different characteristic paths X may have different densities when a path does not 
visit the whole configuration space R^^. The density from the Schrodinger equation is therefore important 
to weight different paths. 

2.4.1. The p error without electrons. Comparing (j2.24[) and (j2.25l) . the difference p — p has contributions 
both from G — G and from the error term M~^/^^^ '^{jp ■ A^aV')- In this section we show heuristically 
how the characteristics can be used to estimate the difference p — p, leading to 0{M~^) accurate Ehrenfest 
approximations of Schrodinger observables 

9{X)p^dX^ j g{X)p{X)dX + 0{M-^), 



in the well known case of no electrons present; Section 17.41 extends this derivation to the case including 
electrons. 

In the special case of no electrons, the X dynamics does not depend on il> and therefore X — X and 
consequently G ^ G. The difference — can be understood from iterative approximations of (I2.13P 

(2.27) _ (t/- Vo)^„+i = ^G^ Ax. (G-iV'«) 

j 

with ipQ = 0. Then = ^ is the Ehrenfest approximation and formally we have the iterations approaching 
the full Schrodinger solution i/)„ — ^ as n — >■ oo. 

In the special case of no electrons, there holds V — Vq, so the transport equation ii/'i = has constant 
solutions; let "01 = 1- Then ■02 — ipi is imaginary with its absolute value bounded by 0{M^^/^); write the 
iterations of ^pn by integrating (j2.27l) as the linear mapping 

n 
fe=0 

which formally shows that 

IV^I^ = + 23fi((0 - V-i) . + - 0i|2 = 1 + O(M-i). 
Consequently this special Ehrenfest density satisfies 

(2.28) p^G-^ ^G'''^tp-i>+0{M~^), 

=p 

since G = G and X do not depend on ip. In the general case with electrons, we show in Section [7]4] that there 
is a solution -0 which is ©(M"^/^) close in L^(dx) to an electron eigenfunction ipo, satisfying 

V{X,-)MX,-) = \oiX)MXr) 

for an eigenvalue Xo{X) £ R and (fixed) nuclei position X. The state ipi equal to a constant, in the case of 
no electrons, corresponds to the electron eigenfunction ipo in the case with electrons present. In the general 
case with electrons. Section [73] shows that still 

(2.29) ijj-iP ^ia + 0{M'^) 

where a — 0{M~^^^) is real and parallel to -00, which implies that the Hamiltonians Hs and He are are 
©(Af"^) close so also G — G — 0{M^^) and consequently the density bound p.28p holds. To obtain the 
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estimate (|2.29p the important new property is that osciUatory cancehation is used in directions orthogonal to 
ipo, reducing the error terms from 0(M~^/^) to 0{M~^), in the case when a spectral gap condition holds. 

2.5. Construction of the solution operator. The WKB-forms p.ip and p.lSp are meaningful when ip 
and ip do not include the full small scale and we verify in Section [7] that both dxip and dx^ are bounded 
independent of AI, in the case of a spectral gap. Section [7] also presents conditions so that ip is 0(M~^/^) 
close to an eigenvector of V in L'^{dx). To replace ip by such an electron eigenstate is called the Born- 
Oppenheimer approximation, which has been studied for the time-independent |50[ [B] and the time-dependent 
[26l |44 l Schrodinger equations by different methods. 

To construct the solution operator it is convenient to include a non interacting particle, i.e. a particle 
without charge, in the system and assume that this particle moves with constant high speed dXl /dt = p\ ^ 1 
(or equivalently with speed one and larger mass); such a non interacting particle does clearly not effect the 
other particles. The additional new coordinate Xl is helpful in order to simply relate the time-coordinate t 
and X^. To not change the original problem (jl.ip . add the corresponding kinetic energy (p})^/2 to E and 
write equation (|2.10|) in the fast time scale r — M^^^t 

3 

and change to the coordinates 

(t,Xo) := {t,XI,XI,X'^,...,X^) e [0,oo) x / instead of (X^ , X^ , . . . , X^) e R^^, 
where X^ = (X^,X^,X|) e R^, to obtain 

(2 30) ^^+2ik^^ =(^-K,)V;-53lGE,A^.(G-i^) 

using the notation w = dw/dr in this section; note also that G is independent of X^. We see that the operator 

2M ^ ^" 



=V-Vo 

is symmetric on L^(R'^'^+'^^~^). The solution of the eikonal equation (|2.4p . by the characteristics 
becomes well defined in a domain U ~ [0,AI^/^t] x R'^^~^, in the new coordinates. Assume now the data 
{Xo,po,za) for Xq G R^^^^ is (iZ)3^-i-periodic, then also {Xr,Pr,Zr) is (LZ)3^"i-periodic. To simplify 
the notation for such periodic functions, define the periodic circle 

T R/(iZ). 

We seek a solution $ of which is (LZ)'^('-'+^'^-'^-periodic in the (x, Xo)-variable. The Schrodinger operator 
Vr has, for each r, real eigenvalues {Xm{T)} with a complete set of eigenvectors {Pm{x, Xq)^} orthogonal in 
the space of x-anti-symmetric functions in L^(T^'^+''^^^^), see [3]; its proof uses that the operator % + -fl 
generates a compact solution operator in the Hilbert space of x-anti-symmetric functions in L^(T^'^+'^^^^), 
for the constant 7 € (0, 00) chosen sufficiently large. The discrete spectrum and the compactness comes from 
Fredholm theory for compact operators and that the bilinear form /rp3(j+jv)-i vVrW + jvw dxdXo is continuous 
and coercive on _ffi(T^('^+^^~^), see [18] and Section [9l We see that V has the same eigenvalues {Am(T)} 
and the eigenvectors {GrPmiT)}, orthogonal in the weighted L^-scalar product 

V ■ w G^^dXo. 

The construction and analysis of the solution operator continues in Sections [7] and [9] based on the spectrum. 
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Remark 2.1 (Boundary conditions). The eigenvalue problem makes sense not only in the periodic 

setting but also with alternative boundary conditions from interaction with an external environment, e.g. for 
scattering problems. The inflow, with data given from the time-independent Schrodinger problem, and the 
outflow of characteristics gives a different perspective on molecular dynamics simulations and the possible 
initial data for the time-dependent Schrodinger equation. 

3. The time-dependent Schrodinger equation 
The corresponding time-dependent Ansatz 

in the time-dependent Schrodinger equation [IH] 

leads analogously to the equations 

dtp + dx^ ipdx^O) = M-^/^Y, • Ax.V)> 
coupled to the Schrodinger equation along its characteristics 

X 



(3.2) 

dx{'>P-vi') 

with the same equation for {X,p,ip, p, z,dxp, dxxp) as for the characteristics p.4p and (|2.13p in the time- 
independent formulation. The Ehrenfest dynamics is therefore the same when derived from the time- 
dependent and time-independent Schrodinger equations and the additional coordinate, introduced by a non 
interacting particle in the construction of the time-independent solution in Section 12. 5[ can be interpreted 
as time. A difference is that the time variable is given from the time-dependent Schrodinger equation and 
implies classical velocity, instead of the other way around for the time-independent formulation. 



4. Surface-hopping and multiple states 

In general the eigenvalue E is degenerate with high multiplicity related to that several combinations of 
kinetic nuclei energy and potential energy sum to E 

-f '^*:TirFy^ '^x^'^dxdX + [ ^*V^dxdX = E [ <^>*<PdxdX, 

with different excitations of kinetic nuclei energy and Born-Oppenheimer electronic eigenstates. When several 
such states are excited, it is useful to consider a linear combination of eigenfunctions X]n=i ' ^" =: $, 

where the individual terms solve (|l.ip for the same energy E. We have 

(4.1) = E^ 

and the normalization X]n=i ll^nllia = 1 implies ||$||i2 = 1. Such a solution $ can be interpreted as an exact 
surface-hopping model. The usual surface-hopping models make a somewhat different Ansatz with the x- 
dependence of ipn prescribed from a given orthonormal basis in L(T'^'^) of wave functions of different energy 
and with explicit time-dependence, see |5H 150) . This section extends the Ehrenfest dynamics to multiple 
states. 
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4.1. Surface-hopping and Ehrenfest dynamics for multiple states. The characteristic {Xn,Pn), the 
wave function V-'n, the density p„ and the phase z„(t) — 0„(X„(t)) determine the time- independent wave 
ftinction ipn and the corresponding Ehrenfest approximation 

(4.2) ELiMrnPny^'e^^'"'"'" =: $ 

yields an approximation to $; the density of state n is now a constant multiple, r„ > 0, of the one-state 
density p„ defined in (|2.26p . normalized to /rpatj+jv) PndX — 1 and X]"=i = ^ the case of no caustics, 
the Ehrenfest states (X„, -0„, /5„), n — l,...,n, satisfying (|2.17p and (|2.26p . are asymptotically uncoupled, 
see Section [SJ a caustic point a is where l//5„(a) = 0. In the case of caustics and n colliding characteristics, 
the phases become coupled and it is necessary to have a sum of n WKB terms for $ to approximate the 
eigenfunction $, in (jl.ip . away from caustic points, see [HI |4T]. Consequently in the presence of caustics, 
surface-hopping approximation may improve poor approximation by Ehrenfest dynamics. 

5. Computation of observables 
Assume the goal is to compute a real-valued observable 

$ . AMX 



for a given bounded linear multiplication operator A — A{X) on L^(T'^^) and a solution $ of (j4.1|) . We have 

The integrand is oscillatory for n ^ m, so that critical points of the phase difference give the main contribution: 
the method of stationary phase (cf. [Ill |41] and Section 17. 5p shows that these integrals are small, bounded 
case when the phase difference has non degenerate critical points (or no critical point) 
and the functions Aipn ■ ipm, On are sufficiently smooth. A critical point a satisfies dx^nici) — dxdm{a) — 0, 
which means that the two different paths, generated by On and Om, through a also have the same velocity p 
in this point. Consequently the critical point must be a caustic point since otherwise the paths are the same. 
That the critical point is non degenerate means that the matrix dxxiOn — ^m)(a) is non singular. 

We see that WKB terms, with smooth phases and coefficients avoiding caustics, are asymptotically or- 
thogonal, so that the density of a linear combination of them, separates asymptotically to a sum of densities 
of the individual WKB terms: 



(5.2) / $ . A^dX = V / Aijn-^ndX + 0{M-'^), 

•^T^" ^1- 

in the case of multiple eigenstates, n > 1, and 



<i> A<^dX = / ^Vi -iJidX 
for a single eigenstate. We will study molecular dynamics approximations of a single state 



(5.3) / AiPn-^ndX^ A{X)pn{X)dX. 

in the next section. 

In the presence of a caustic, the WKB terms can be asymptotically non orthogonal, since their coefficients 
and phases typically are not smooth enough to allow the integration by parts to gain powers of Af^^/^. Non 
orthogonal WKB functions tells how the caustic couples the WKB modes. 

To numerically compute the integral (15.31) requires to find an approximation /)„ of p„, using the Ehrenfest 
solution (X„,^„), and to replace the integrals by quadrature (with a finite number of points X). The 
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quadrature approximation is straight forward in theory, although costly in practical computations. Regarding 
the inflow density p„ | ^ there are two situations - either the characteristics return often to the inflow domain 
or not. If they do not return we have a scattering problem and it is reasonable to define the inflow-density 
pn\i as an initial condition. If characteristics return, the dynamics can be used to estimate the return-density 
f}n\j as follows: assume that the following limits exist 

(5.4) YimT^oo ^ !^ A{Xt)dt = J^,^ A{X)pn{X)dX 

which bypasses the need to find pn\j and the quadrature in the number of characteristics. A way to think 
about this limit is to sample the return points Xt (z I and from these samples construct an empirical return- 
density, converging to /5„ | ^ as the number of return iterations tends to infinity. We shall use this perspective 
to view the Eikonal equation (j2.16l) as a hitting problem on /, with hitting times r (i.e. return times). We 
allow the density p\i to depend on the initial position Xq; the more restrictive property to have p\i constant 
as a function of Xq is called ergodicity. 

6. Approximation of a deterministic WKB state by Hamiltonian dynamics 

A numerical computation of an approximation to /x3n • AipndX has the main ingredients: 

(1) to approximate the exact characteristics by Ehrenfest characteristics (|2.17|) . 

(2) to discretize the Ehrenfest characteristics equations, and either 

(3) if p\j is an inflow-density, to introduce quadrature in the number of characteristics, or 

(4) if is a return-density, to replace the ensemble average by a time average using the property (|5.4p . 
This section presents a derivation of the approximation error in step one, in the case of a return density, 
avoiding the second, third and fourth discretization steps studied for instance in pi [TUl 

6.1. The Ehrenfest approximation error. This section shows that the Ehrenfest system, written as a 
Hamiltonian system (jl.7p . approximates Schrodinger observables without caustics. We see that the approxi- 
mate wave function $ defined by 

(6.1) ^ ^ p'/^e'^'"''\ 

where ij) = ^e**^''"" -^o '* '^''^(«)''« and (X,0) solves the Hamiltonian system ([2:2T|) or the system ([2:20|) . is an 
approximate solution to the Schrodinger equation (jl.ip 

(6.2) {H - E)^ ^ E ^P""^)^ 

3 

since by ^T^, ([^T^ . ^^IJ^ and (PTf|) 

{H -E)^ =- {iM-^/"^^ -{V- Vo)^) /5i/2e'M 



l/2fi 



(6.3) 



+ +^.V^-E) pV2^e'MV2fl 



Therefore $ approximates a single WKB eigenstate $, satisfying iJ$ = The following theorem presents 
conditions for accurate approximation error of observables 



/ g{X) $ • $ - / g{X) $ • $ dX, 

J-pSN Jt^N '^—^ 



=p(X) =p{X) 
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using that either the spectral gap condition (|6.4|) or the crossing eigenvalue condition ()6.5p holds, based on 
the electron eigenvalues A„, defined by V{X)^pn{X) — A„(X)?/;„(X) in L?{dx). 

The spectral gap condition. The electron eigenvalues {A„} satisfy for some positive c the spectral gap 
condition 

(6.4) inf |A„(y) - Ao(r)| > c, 

n^Q, YeD 

where D :— {Xt : t > 0}U {Xt : t > 0} is the set of exact and approximate nuclei positions. 

The crossing eigenvalue condition. Assume that there is a positive constant c such that the sets of eigenvalue 
number n and eigenvalue crossing times a, with maximal hitting times r (defined in ()7.11|) ). 

, . nx ^{{n,a) : A„(X,) = Ao(X,), 0<CT<T,n^O} 

% ={(n,a) : A„(X,) = Ao(X,), 0<a<T,n^O} 

consists of finite many crossing times a with precisely one n value for each crossing time (but possibly different 
n at different crossing times) and 



|^(A„(X,)-Ao(X,))|(„.*)e^,, >c 
|^(A„(X0-Ao(lt))lM)e^2^ >c. 

Theorem 6.1. Assume that the electron eigenvalues have a spectral gap (j6.4p and a = I — S, or that they 
form non degenerate critical points (|6.5|) and a = 3/4, then Ehrenfest dynamics (I2.23|) . assumed to avoid 
caustics and have bounded hitting times r in (|7.1ip . approximates time-independent Schrddinger observables 
with the error bounded by 0{M~"') : 

(6.6) /t3« g{X)p{X)dX = /T3iv g{X)p{X)dX + 0{M-"), for any S > 0. 

The proof is in Section [71 with the main steps separated into six subsections: spectral decomposi- 
tion, Schrodinger eigenvalues as a Hamiltonian system, stability from perturbed Hamiltonians, the Born- 
Oppenheimer approximation, the stationary phase method, and error estimates of the densities. The ob- 
servable may include the time variable, through the position of a non interacting (non-charged) nucleus 
moving with given velocity, so that transport properties as e.g. the diffusion and viscosity of a liquid may by 
determined, cf. [13]. A case with no caustics is e.g. when mmx{E — Vo{X)) > in one dimension X eR. 

6.2. The Born-Oppenheimer approximation error. The Born-Oppcnhcimer approximation leads to the 
standard formulation of ab initio molecular dynamics, in the micro-canonical ensemble with constant number 
of particles, volume and energy, for the nuclei positions X, 

(6.7) Xt =pt 

by using that the electrons are in the eigenstate ipo with eigenvalue Aq to V, in L^{dx) for fixed X. The corre- 
sponding Hamiltonian is Hbo{^ ,P) '■= |pP/2-l- Ao(A"). The following theorem shows that Born-Oppenheimer 
dynamics approximates Schrodinger observables, based on a single WKB eigenstate, as accurate as the Ehren- 
fest dynamics, using the approximate Schrodinger solution 

with density p := $ • $ and G defined by dlogG/dt — divp/2, as for G in (|2.9p 

Theorem 6.2. Assume that the electron eigenvalues have a spectral gap (j6.4l) and a = 1, or that they form 
non degenerate critical points (|6.5p and a = 1/2, then the zero-order Born-Oppenheimer dynamics (j6.7p . 
assumed to avoid caustics and bounded hitting times r in (j7.1ip , approximates time-independent Schrddinger 
observables with error bounded by 0{M^"^^) : 

(6.8) /t3« giX)p{X)dX = /^3« g{X)p{X)dX + 0{M-"+'), for any 5 > 0. 
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The proof is in Section [71 

Remark 6.3 (Why do symplectic numerical simulations of molecular dynamics work?). The derivation of 
the approximation error for the Ehrenfest and Born-Oppenheimer dynamics, in Theorems 16.11 and 16.21 also 
allows to study perturbed systems. For instance, the perturbed Born-Oppenheimer dynamics 

p =-dxXo{X)-dxH\X,p,^j), 
generated from a perturbed Hamiltonian Hbo{XtP, "0) + H^{X,p, tp) — E, with the perturbation satisfying 
(6.9) \\H^\\l-<'<S for some 5 e (0,cx3) 

yields through (|7.13p and (|7.16p an additional error term 0{5) to the approximation of observables in (j6.6p . 
So called symplectic numerical methods are precisely those that can be written as perturbed Hamiltonian 
systems, see [46j . and consequently we have a method to precisely analyze their numerical error by combining 
an explicit construction of with the stability condition (|6.9I) to obtain 0{M~^+S) accurate approximations. 
The popular Stormer-Verlet method is symplectic and the positions X coincides with those of the symplectic 
Euler method, for which is explicitly constructed in |46) with S proportional to the time step. 

7. Analysis of the molecular dynamics approximation 

This section continues the construction of the solution operator started in Section [2751 It is written for the 
Schrodinger characteristics, but it can be directly applied to the Ehrenfest case by replacing V hy V — Vq (by 
formally taking the limit Af — > oo in the term {2M)^^GJ2j ^Xj{G~^ip))- Assume for a moment that V is 
independent of r. Then the solution to (|2.30p can be written as a linear combination of the two exponentials 

Ae'^°'+ + Be"°'- 

where the two characteristic roots are 

a± = (p})^(-l±(l-2(p})-^F)V2). 
We see that e*'^"- is a highly oscillatory solution on the fast r-scale with 

a^ = -2{pl)^ + V + 0{vy{pir), 

while 

(7.1) a+ = -V + 0{Vy{pir). 
Therefore we chose initial data 

(7.2) #|r=o = -a+iJAr^a 

to have B — 0, which eliminates the fast scale, and the limit — > oo determines the solution by the 
Schrodinger equation 

iijj = Vip. 

In the case of the Ehrenfest dynamics, this equation with V replaced by — Vq is the starting point. The 
next section presents an analogous construction for the slowly, in r, varying operator V. 



-1 
2{p\fV -2{p\f 
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7.1. Spectral decomposition. Write (|2.30|) as the first order system 

iij} = V 

V =2{p\fi{V^-v) 

which for -0 := {ij:, v) tafces the form 

= iA'il) A := 

where the eigenvalues A± , right eigenvectors q± and left eigenvectors q^"^ of the real "matrix" operator A are 

q+ 
q- 

We see that A+ = ~V + and A_ = -2{p\f + V + 0{{p\y^). The important property here is 

that the left eigenvector limit \m\pi^rx, = (1, 0) is constant, independent of r, which implies that the q^ 
component q'^^'ip — ip decouples: we obtain in the limit — )■ oo the time-dependent Schrodinger equation 

= -q^ Ar^pr _ 

= -A+(r)Vi(r) 
= Vr^'ir), 

where the operator Vr depends on r and {x,Xq), and we define the solution operator S 

(7.3) ^(r) = ^.,0^(0). 

As in (|7.2p we can view this as choosing special initial data for "0(0); from now on we only consider such data. 
The operator V can be symmetrized 

(7.4) Vr := G;'VrGr = {V- Vo)r - ^jly E ^xi^ 




J 

'2( 



with real eigenvalues {Am} and orthonormal eigenvectors {pm} in L" [dxdXf)) , satisfying 

VrPmir) = Am(r)p™(T), 

see Section ini Therefore Vr has the same eigenvalues and the eigenvectors Pm := GrPm , which establishes 
the spectral representation 

(7.5) 14^(.,r,.) = VA„(t) / ^{■,0,-)-p^Gr^dXoP^{T), 
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where the scalar product is 

(7.6) [ i^-p^G-^dXo. 

We note that the weight G^^ on the co-dimension one surface T^^~^ appears precisely because the operator 
V is symmetrized by G and the weight G~^ corresponds to the Eulerian-Lagrangian change of coordinates 
(|2.12p . The existence of the orthonormal set of eigenvectors and real eigenvalues makes the operator V 
essentially self-adjoint in the Lagrangian coordinates and hence the solution operator S becomes unitary in 
the Lagrangian coordinates. In the case of the Ehrenfest dynamics the weight is the density p = G~^; 

7.2. Schrodinger eigenvalues as a Hamiltonian system. Our error analysis is based on comparing 
perturbations of Hamiltonian systems. This section establishes the Schrodinger eigenvalue characteristics as 
a Hamiltonian system, analogously to the Hamiltonian system (|2.23l) for Ehrenfest dynamics. 

To write the Schrodinger eigenvalue characteristics as a Hamiltonian system requires to make an asymp- 
totically negligible change in the definition of Vq in (j2.10p : the Eikonal equation used the leading order term 
Vq = ip ' Vip/^ ■ tjj. An alternative choice is to replace V by the complete Schrodinger operator 

Vi> := yVi-G(2M)-i^AjfjV'G~^) 

n 

and symmetrize 

(7.7) vo:= ^-n+n-^ 



2 J^P-tPG-^dX/ jG-^dX 

The integral in X is used in the denominator to avoid the non constant ip ■ ip; since V is essentially self-adjoint 
in the weighted space L^{G~'^dX) in (|7.6p . corresponding to Lagrangian coordinates, the weighted norm 
in the denominator is constant in time and consequently the denominator can be treated as a constant when 
differentiating the Hamiltonian to obtain the Hamiltonian system. 

This definition of Vq has the same leading order term as the previous (|2.19|) . Equation ()2.10|) showed that 

={H-E)^ 

= (( - + - ^o)Vi - 217 GAx. (^G-i))G-i 

for any definition of the scalar Vqj which we now use for the new Vq. Note that, it is only when {V — Vo)ip 
is small that we expect the X-derivative of 'ip to be bounded, since otherwise the e**^ ^ scale will pollute 
derivatives: p ■ dxip = —iM-^/^{V — Vo)'ip. Therefore, the interesting case is when -0 is close to an electron 
eigenfunction and Vq is near the corresponding eigenvalue. 

The goal here is to show that the Eikonal equation, corresponding to the second term above, becomes a 
Hamiltonian system for the value function 6{X^p,ip), so that it also generates the Schrodinger equation in 
the first term above, analogously to the Ehrenfest dynamics (|2.23p . Define (p = 2M^^'^(p and the Hamiltonian 

(7-8) Hs{X,p,ip) -.^^ — + , 

similar as in (j2.23p . which yields the Hamiltonian system 

X ^p 



(7.9) p 



4>-dxV<i>+dxV4>-4> 
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using that V is symmetric in the weighted scalar product L^{dxG^^dX), see (|7.4p . The phase factor change 



P 2 / ijj-tp G-^dX/ f G-^dX 

i^P = Ml/2 _ ^ 

for 1^0 defined by (fTT)) and 

The last two equations imply that {X,p,^) and the weight G, defined by (|2.9p (using the new Vo)j generate 
a WKB function $ := G~^'>pe^^^^^^ that solves the Schrodinger equation = in (jl.ip . Therefore, the 
Hamiltonian (|7.8p generates an exact solution to the Schrodinger equation. 

7.3. Stability from perturbed Hamiltonians. This section derives error estimates of the weight functions 
G when the corresponding Hamiltonian system is perturbed. 

To derive the stability estimate we consider the Hamilton- Jacobi equation 

in an optimal control perspective, with the corresponding Hamiltonian system 

Yt ^dqH{quYt) 
qt ^-dYH{qt,Yt). 

Define the "value" function 



9{Yo) = 9{Yt)- / hiq,,Y,)ds 
Jo 

where the "cost" function defined by 

:^q>dgHiq,Y)-H{q,Y) 
satisfies the Pontryagin principle (related to the Legendre transform) 
(7.10) H{q, Y) = sup {q . dgHiQ, Y) - hiQ, Y)) . 



Let 



be defined by the hitting problem 



e{Yo)^e{Y,)- / %„n)ds 

Jo 

using the hitting time r on the return surface / 

(7.11) T inf{t : Yt e I kt>0}. 

Define analogously for a perturbed Hamiltonian H, the dynamics {Yt,qt), the value function 9 and the cost 
function h. 

We can think of the difference 9 — 9 as composed by perturbation of the boundary data (on the return 
surface /) and perturbations of the Hamiltonians. The difference of the value functions due to the perturbed 
Hamiltonian satisfy the stability estimate 

(J ^2) 9{Yo)~9{Yo) >9{Yr)-9{%)+^^{H~H){dY9{Yt),Yt)dt 

^' ' 9{Y^)-~9{Yo) <9{Yr)-9{Yr)+J^{H~H){dY9{Yt),Yt)dt 
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with a difference of tlie Hamiitonians evaiuated along tlie same solution path. This result follows by dif- 
ferentiating the value function along a path and using the Hamilton- Jacobi equations, see Remark 17.11 and 

Assume that 

(7.13) sup \iH - H){q,Y)\ ^ 0{M-"). 

The stability estimate, for all characteristic paths and following subsequent hitting points on / (as in the 
discussion after (|2.6p ). then yields the same estimate of the difference in the boundary data 

provided the maximal hitting time r is bounded, which we assume. If the return surface is the plane for 
which the Xi particle has its first component equal to its initial value, the hitting time is the time it takes 
until this particle has the first component equal to that initial value again; since the X-dynamics does not 
explicitly depend on M it seems reasonable that one can find a return surface such that the hitting times are 
bounded. Next, the representation can be applied to interior points to obtain 

\\9-9\\l^ ^0{M-''). 

When the value functions 6 and 6 are smoothly differentiable in X {X is a part of the Y coordinate) 
with derivatives bounded uniformly in M, the stability estimate (|7.12l) implies that also the difference of the 
second derivatives has the bound 

||^Ax„6'-^Ax„6'||l- =0(Af-"+*), for any 5 > 0. 

n n 

We will also use that 

(7.14) M-i/2!3(V' • J2 ^x,A)/{'^ ■ i') = C(A/"") 

n 

which is verified in Section 17.4.41 

Our goal is to analyze the density function p satisfying the convection equation (j2.9p . i.e. 

(7.15) dxe •dx\og p^~^l:^xj + d 

n 

where the function d — dso^ d — dE ot d — ds with 

dso d_E = 
ds - -M-i/2cj(7A • ^xA)/{i> ■ 4') 

n 

in the case of Born-Oppenheimer, Ehrenfest or Schrodinger dynamics, respectively. We have dxO»dx^ogp — 
d\og p{Xt) / dt and we interpret equation (|7.15p as a Hamilton- Jacobi equation and similarly for log p. Then 
the stability of Hamilton- Jacobi equations, as in (j7.12p . can be applied to (j7.15p . with the Hamiltonian 

Hp{dx \ogp,X) dx0{X) . dx logp(X) + ^ AxJ{X) ~ d, 

n 

to obtain the sought error estimate 

(7.16) \\\og p- log p\\l^^O{M-^+'). 

In this sense we will use that an 0{M~°') perturbation of the Hamiltonian yields an error estimate of almost 
the same order for the difference of the corresponding densities p — p. 
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The Hamiltonians we use are 

Hs ^^■£^^+\{4>-V{X)4> + V{X)4>-^)-E 

He =^-^ + (t>-V{X)(t>- E 

Hbo -^ + Ao(X)-£; 

with the cost functions 

hs - + ^ + i (<^, • V{X)^, + V{X)4>, ■ ^, - ■ V{X)cPr - V{X)^r ' ^r) 
He = S + ^ + • - <Pr ■ V{X)(t>r) 

hso =S+^-Ao(X) 

and the primal and dual variable {Y]q) = (X, 2^/^M~^/^(/)r;p, 2^/^M~^/'*(/)i) in the case of Schrodinger and 
Ehrenfest dynamics, as in (|2.22l) . and the variable {Y;q) — {X,{tpQ)r]p,{^o)i) for the Born-Oppenheimer 
dynamics. For the Born-Oppenheimer case the electron wave function is the eigenstate V'o! one can identify 
the Born-Oppenheimer dynamics with an Ehrenfest dynamics where the mass is set to infinity, since in this 
limit the Ehrenfest wave function becomes the eigenfunction, see Section [7.41 

Remark 7.1. This remark derives the stability estimate (I7.12p . The definition of the value functions imply 

e{%)~ I h{Qt,Yt)dt-{diYr)- I h{quYt)At) 



- - /; h{qt,Yt) dt + 6{Yf) - e{Ya) - e{%) 

(7.17) e(%) 

= - /; h{quYt) dt + /; de{Yt) + e{%) - e{Yr) 

= /; -HquYt) + dyOiYt) . d^H{qu Yt) dt + e{Yf) - 0{%) 

" V ' 

<H{dYe{Yt),Yt) 

< j^iH - H){dYe{Yt),Yt) dt + e{Yr) - e{Yr), 

where the Pontryagin principle (|7.10l) yields the inequality and we use the Hamilton- Jacobi equation 

H{dYe{Yt),Yt)^0. 

To establish the lower bound, replace along Yt by 6 along Yt and repeat the derivation above. 
7.4. The Born-Oppenheimer approximation. 

7.4.1. With an electron eigenvalue gap. To better understand the evolution (|2.17p of ■0 and estimate the error 
term in (|6.3p . we use the decomposition ■(/' = fAo © where "00(0 an electron eigenvector of Vt = V{Xt), 
satisfying Vt^o{t) — Xo{t)'ipo{t) in L?{dx) for an eigenvalue \o{t) g R, and is chosen orthogonal to tpo in 
the sense • V'o = 0. We assume that the electron eigenfunction ijjviiX) : T'^'^ — > R is smooth as a function 
of X. This Ansatz is motivated by the residual 

(7.18) i?^o V-o + iM^/^{V - Fo)V'o 

being small. Below we verify that if)-^ is 0(Af ~^/^) with a spectral gap assumption and 0{M~^/'^) for a case 
without a spectral gap; in the case of a spectral gap, this estimate yields 
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Since 



for = 0{1), and 



(ipo + V' ) • wo + V' ) 



The Schrodinger equation Rip = in (|2.17p implies that the perturbation satisfies the following 
Schrodinger equation with the source term iRipQ 



(7.19) 

and the solution representation 
(7.20) 



for the solution operator ipt =: Stfii'o- Split the source term into its projection on -00 and its orthogonal part 

i?'0o • -00 



00 ® RlPo 



and note that it is enough to study the projected source term RipQ to determine ip-^ on the orthogonal 
complement of ipQ . Integrate by parts to obtain 



(7.21) 



= St,s 



StAV~Vo)-^RiP^{.s) 



If there is a spectral gap 
(7.22) 

for a positive constant c, we have 



inf \K{Xt)-\o{Xt)\>c, 



{V - Voy^RlP^ - - ^0)"Vn R^O ■ V'n = I] (An ^ K))"Vn i?0'o^ ' 0n = ^(1) 

n>0 nT^O 

We have by (j?:^ and (j?:^ 

^-^(0 = ^t,o(0^(O) - iA/-i/2(F - Vo)-'RiP^{0)) + iM-^/\V - Vo)-'RiP^{t) 

- r 5t^,^((F - Vo)-'RiP^){s)ds, 



where the last integral can be integrated by parts again as in (|7.2ip to reduce the integral with a factor 
By choosing V'-^(O) = iM~^/^{V - Vn)''^ Rip^ (0) we have 

(7.23) iP^{t)=iM-'^\V-Vo)-'R'iJj^{t)+iM-'/^ j" St^^±(^v - Vo)'' Rii4){s)ds . 



0{M- 



MOLECULAR DYNAMICS DERIVED FROM THE SCHR.ODINGER EQUATION 



25 



We may use that the paths Xt return to the co-dimension one "return" surface /, defined in p.6|) . after time 
i = which by (fT^ then yields 

(7.24) II / St,sRi^^{s)ds\\L2 - 0(M-i/2), 

Jo 

and if longer return times are needed the integration by parts can be iterated k times to gain a factor of 
TVf in the integral. The function ip-^ is then determined by (|7.20p from the initial X € I and its successive 
returns to the set /, similarly as the phase function 9 in (|2.6|) . Therefore, the unitarity of St,Q together 
with the bound (|7.24p imply in (|7.20p that the Ehrenfest approximation = ipo + w satisfies 

(7.25) II^^IU. =0(A/-i/2), 

where (|7.23|) shows that the error in different hitting time intervals do not substantially accumulate in time. 
Since the function f/'o depends smoothly on X and a derivative of the integral in ()7.23p yields at most a factor 
Af^/^ (when differentiating the solution operator 5), equation (j7.23p shows that the derivative with respect 
to the Lagrange coordinates satisfies 

lia^^lU. = o(M-V2), 

and integrating by parts once more in (|7.23p also gives 

(7.26) ||a2^/.^||i2 =0(A/-i/2). 



7.4.2. With crossing electron eigenvalues. Since the solution operator is unitary it can be written as a highly 
oscillatory exponential: St^s — e**^^^^'3'^'*\ where Q is an essentially self-adjoint operator, cf. [3]. Therefore, 
the method of stationary phase, cf. [41 and Section [731 can be applied to extend (|7.2ip to the case when 
~ has finitely many non degenerate critical points (7^, where the eigenvalue crossings satisfy 

(7.27) \'^{^n{s) — Ao(s))|s=crfc > c for a positive constant c and n^Q. 

With this assumption, the method of stationary phase, see Section [7.51 implies 



f St,sRi^^{s)ds = 0{Ar 
Jo 



1/4n 



and 

(7.28) IIV-^ll =0(Af-i/4), 

since by letting the hitting surface / include the points where the electron eigenvalue cross we have 

^^(t) = 5t,o(^^(0) + Af-i/V) 

for a smooth and bounded function r oi X and for times t close to the time of a path hitting crossing 
eigenvalues. To obtain bounds on the derivatives, we note that if the initial data on the hitting surface is 
chosen as ^"^(0) = -M~^/'^r + 0{M-^/^), we obtain 

(7.29) liaVlU^ -o(M-'/"), 

since two derivatives of the solution operator increases the result by at most a factor 0{M). 
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7.4.3. Difference of Hamiltonians. To estimate the difference in the Hamiltonians we use the bounds (|7.26p 
and (|7.29p on derivatives 

^•G^Ax,(G-iV')-0(l), 
in the Lagrange coordinates. We conclude that 

The largest bound 0{M-^/^) for crossing eigenvalues can possibly be improved, since when the eigenvalues 
cross we did not use the smallness of ■ {V — Vo)^''" in the neighborhood of crossing eigenvalues but 
only the size of ^jJ'^ = 0{M~^/'^): in a neighborhood of the crossing point the function ip-^ increases from 
0{M~^/'^) outside to C(M~^/^) inside, for the component corresponding to the crossing eigenvalue, while 
{V — Ao)V'^ vanishes asymptotically at the crossing; a careful balance of the size of this neighborhood and of 
{V — \o)^^RiIjq outside could give a better estimate of 'tp'^ • [V — ^o)'4'^ ■ 

7 A A. Estimates of the Hamiltonian for the density. We have ■0 = G^^ii'o + V-"^)- To conclude that (|7.14p 
holds, we split into three terms where the first is 

9(G-Vo-^Ax„(G-Vo)) =0 

n 

since all functions here are real valued; then we have 

(7.31) M-/.c,(G-.^„.^^^^(G-.^.) + C-.^,..^^^^(C-.^„) r O^iMp s,^^,., 
and 

^ ^ 0(M ^) crossmgs. 

by the derivative bounds (|7.26p and ()7.29p . We conclude that 

7.5. The stationary phase method. 

7.5.1. Real valued phases. The theory of semi-classical approximations [411 114) use the method of stationary 
phase to reduce the study of oscillatory integrals 



ts 



(7.33) / e^^'^''''3(^'/(s)d, 

JR 

to points a where the real valued phase function Q is stationary, i.e. Q'{(j) = 0, as follows; such integrals 
over domains away from the critical points Ufc, yields by the integration by parts (I7.2ip contributions of size 
0{M~-^/^). In a neighborhood of a critical point a = 0, the phase can by Taylors formula by written 

Q{s) = Q(0) + - / Q"{ts)dt s^. 
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Introduce the change of variables s :— s{Q{s)/Q{0)y^^ and its inverse map s{s) to write 

/>oo 

ga/V=,= Q(,)/2 ^^^^^^^ ^ / ^M^/^S^QiO)/2 f^s{S))s'{S) ds 



The corresponding error term based on the integral over the domain (— oo,0) defines / on (— cx),0). Fourier 
transformation J-, following |14| , shows the expansion 

J gjMi/^s2g"(o) /(s)rfs 

^jR^{e'^'''''''^"^°^){uj)T-^ (/(s)) {Lj)duj 
(7-34) jr^(2^)i/2M-i/4|g"(o)|-i/2e»-gnQ"(o)/4g-,:(Q")-iAf-V2„V2j-i(/(s))(^)rf^ 



(27r)V2Af-V4|Q"(0)|-l/2e<-gnQ''(0)/4^-^M_j:^(^g,,(0)-1^2)fcj(^^ 



The original integral (j7.33l) can by a smooth partition of unity be split into a sum of the integrals over domains 
containing only one or no critical point. 

The stationary phase method can be extended to oscillatory integrals in R'^^, with the factor AI^-^^'^ 
replaced by M~^'^/^ and |g"(0)| replaced by | detg"(0)| in (TTTMl) . cf. [Hill]. 

7.5.2. Self-adjoint phases. In this section we extend the stationary phase method to the phase generated by 
the solution operator S. Since S is unitary, we can write the solution operator St.s = e'*^ ' Qt,s ^ where Qt,s is 
essentially self-adjoint in L^{dx). Therefore Qt^s has real eigenvalues {A„(s)} with corresponding orthonormal 
real eigenfunctions {pn{s)} (satisfying Qspn{s) = A„(s)p„(s)) that transforms the operator valued phase to 
a sum of real valued phases 

/* St,sRsds = /* e^'^^'^Q^Rsds = y f e'Af^''^"(^) (p„(s) • i?,) p„(s) ds 

Jo Jo n J° " « ' 

= (p„(t)-fl,) p„(t) 

and the stationary phase method can now be applied to each term separately. 

Note that pn{t) does not depend on s, since the data for p„ is that they are the eigenfunctions of Qt^t = 
F — Vb at time t: the eigenvalue relation 

and {t — s)^^Qt,s V — Vq as s —i' t shows that Pn{t) are the eigenfunctions to 1^ — Vq; let /„(t) :— 
e**^^^^^"^'*^p„(s), then the continuity at t implies = Pnit), so that Pnit) — e**'-^^^^^"('*'*)p„(s; i), which 

yields (p„(s) • i?^) p„(s) = {p„it) ■ Rs) p„{t). 
The definition of Qs implies that 

fQs ^V{Xs)-Vo{Xs) 

hQs ^Uv{x,)-v,{x,)), 

and by differentiation of the eigenvalue relation we obtain 

-7-A„(s) = pn ■ {V - Vb)p„(s). 
ds 

When i — >■ s we obtain A„ ~ A„ and A„ = A„, so that the condition (|6.5I) for crossing eigenvalues of A„ yields 
the condition for non degenerate critical points of A„ . 
Our analysis shows that 

Po{t) = e^^^'''^"(^^*)po(s) = St,sPo{s) 
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which by (|7.28|) satisfies po — i^o + 0{M ^/'^). Analogously, with a similar assumption on non degenerate 
critical points in A„i — A„, for all m and n, we obtain for the other eigenfunctions e*^^^ ^ -f^ ^"^^''~^°^^'''^^pn{s) — 
ipn + ©(A/^^/^), and we can conclude that p„ is 0{M'^/^) close to the electron eigenfunction 



7.5.3. Analysis of the Schrodinger dynamics. The analysis above can be applied to the exact solution operator 
St^s — 6**^^^^*^= and V replacing the Ehrenfest operator 5* and V — Vq. Let ip = ipo ® i''^^ then -ijjt = Stfiipo 
and R^po = "00 ^ M^I'^V-ipQ. Integrate by parts to obtain 



\l St.sRi=k{s)ds - ll ±(St,s)tM~'/^V-' R^^{s)ds 



(7.35) 



and use 



= St.s 

St,sV-^R^^is)' 



-^^'^y^jls,,,±{v-^Ri^^){s)d. 



s 



V-^RiP^ ={l+{V-V^Vo){V-Vo)-^) \V-Vo)-'Ri^^ 

= En^O {I+{V-V~ VoKV - Vo)-'y\Xn - Vo)-'iJnR^4 ' 0™ 

to deduce as before that -0^ = ©(M"^/^) when the spectral gap condition holds and V'n are smooth. 

In the case of crossing eigenvalues, the method of stationary phase, projected to the eigenfunction of 
the singularity, is applicable with S, in the weighted scalar product / v ■ wG~'^dX where the symmetric 
generator is {V — Va)G~^ip — (2M)~^ J2n {G~^ip)- In the new variable ip — G~^ip the operator becomes 
V — Va — {2M)^^ J2n with its time derivative equal to (i(V^—Vo) /(it, as in the Ehrenfest case. The method 
of stationary phase for 5 can therefore be used to study the behavior of the oscillatory integral /J St^sRipo 
in the subspace of a single eigenfunction ipn, in a neighborhood a critical point where A„ — Vq vanish, n 7^ 0; 
in the orthogonal subspace to ipn and ipo, the oscillatory integral can be estimated by integration by parts as 
in (17.351) . With such a projection into a single eigenfunction direction, the source term is i?V'o" ' V"" and the 
non degeneracy condition for the stationary phase method becomes 

■ Q{0)4>n\ = • d{V - Vo)/dt4,r,\ = |A„ - T>o | > C. 

Note that the projection to a one dimensional subspace is possible when the critical points of the different 
eigenvalues are separated. If two or more eigenvalues have the same critical point, a more careful study would 
be needed. 

In conclusion we have the same estimates for ip —: ipo Q) ip^ as for ^p: 

(7.36) IIV'^IU^ = 0(Af-i/2) 
with a spectral gap, respectively 

(7.37) |1V>||l2 =0(Af-i/4) 

with non degenerate eigenvalue crossings, and as in (j7.3ip . we also obtain 



(7.38) 



1 r ^ CG^^i/j) — / ^^^'^ with a spectral gap 

2M ^ ^ 1 OiM^^) with eigenvalue crossings. 
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7.6. Error estimates for the densities. The densities are 



ps = G^'^ip ■ ij} for the Schrodinger equation 

Pe ~ G^^ for the Ehrenfest dynamics 

Pbo = for the Born-Oppenheimer dynamics. 



We have from (TflSOl) . (17.36117.381) and (061) 

G-^ = G-2 +C'(Af-") 
= G-2 +C'(Af-") 



and by (|7.361l7.37p 

(7.39) ^■^ = 1 + 0{M~"), 




The sharper estimate for the density obtained using (j7.16p and (|7.32p . instead of (|7.39l) . improves the error 
bound to 0(M-3/4) ([7:39|) foj. the case of Ehrenfest dynamics approximating Schrodinger observables with 
crossing electron eigenvalues. 



In this section we analyze a situation when the WKB eigenstate, (corresponding to a degenerate eigenvalue 
of the Schrodinger equation (jl.ip ) is perturbed from the corresponding electron ground state by thermal 
fluctuations, which will lead to molecular dynamics in the canonical ensemble with constant number of 
particles, volume and temperature. To determine the stochastic data for the Schrodinger wave function 
in (|7.9|) requires some additional assumptions. Inspired by the study of a classical heat bath of harmonic 
oscillators in [52] . we will sample randomly from a probability density given by an equilibrium solution 
f (satisfying dtf = 0) of the Liouville equation dtf + dp^Hsdrs f — drgHsdpg f = 0, to the Schrodinger 
Hamiltonian dynamics (|7.9p . There are many such equilibrium solutions, e.g. / = h(Hs) for any differentiable 
function h and there may also be equilibrium densities that are not functions of the Hamiltonian. The 
Hamiltonians for the Schrodinger and Ehrenfest dynamics differ by the term 



which is small of is smooth enough, as was verified in Theorems 16.11 and 16.21 Since the Ehrenfest Hamiltonian 
is simpler to understand, in particular showing that smooth are most probable, we first focus on the 
Ehrenfest dynamics. 

To find the equilibrium solution to sample from, we may first consider the marginal equilibrium distribu- 
tion for the nuclei in the Ehrenfest Hamiltonian system. The equilibrium distribution of the nuclei is simpler 
to understand than the equilibrium of electrons: in a statistical mechanics view, the conditinal probability of 
finding the now classical nuclei with the energy He ■= 2^^\p\'^+(p-V{X)(p, for given (p, in (|2.22p is proportional 
to the Gibbs-Boltzmann factor exp(— if^/T), where the positive parameter T is the temperature, in units of 
the Boltzmann constant, cf. [ID], [30j- Assuming that the equilibrium density is a function of the Hamiltonian 
and that the marginal distribution for p is the Boltzmann distribution proportional to e^l''' /f^^), its marginal 
distribution therefore satisfies 



8. Approximation of stochastic WKB states by stochastic dynamics 



1 

M 



5R(0-G^Ax„(G-V)) 



n 
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for some function C : R'^^ — >■ R and Fourier transformation with respect to \p\ G R of this equation imphes 

the Gibbs distribution 

KHe) = ce-"-'^, 

for a normahzing constant c = 1/ J exjp(—HE/T)d(j)rd(j)idXdp. In this perspective the nuclei act as the 
heat bath for the electrons. The work [33* considers Hamiltonian systems where the equihbrium densities 
are assumed to be a function of the Hamihonian and shows that the first and second law of (reversible) 
thermodynamics hold (for all Hamiltonians depending on a multidimensional set of parameters) if and only 
if the density is the Gibbs exponential exp(— iJ^/T) (the "if part was formulated already by Gibbs |25]); in 
this sense, the Gibbs distribution is more stable than other equilibrium solutions. An alternative motivation 
of the Gibbs distribution, based on the conclusion of this work (in a somewhat circular argument), is that 
the nuclei can be approximated by classical Langevin dynamics with the unique invariant density exp ( — 
|pp/2 — Ao(X))/T, which is an accurate approximation of the marginal distribution of exp(— iJ^/T) when 
integrating over all the electron states 0, see Lemma lOl and Theorem l8.3l Note that there is only one function 
of the Hamiltonian where the momenta pj are independent and that is when f{rE,PE) is proportional to 
exp(-iJB/r). 

Since the energy is conserved for the Ehrcnfest dynamics -now viewed with the electrons as the primary 
systems coupled to the heat bath of nuclei- the probability of finding the electrons in a certain configuration 
(j) is the same as finding the nuclei in a state with energy He^ which is proportional to exp(— Jf^j/T) in the 
canonical ensemble. This conclusion, that the probability to have an electron wave function is proportional 
to exp{—HE /T)d(p'^d(p^ is our motivation to sample the data for (p from the conditioned density generated 
be exp {—(j) ■ V{X)(f)/T)d(t/'d<l)^: since we seek data for the electrons, we use the probability distribution for 
conditioned on {X,p). 

In (O) 

we esta-blislied the Rsymptotic density Pn^n 

for a multiple state system with probability r„ for 
eigenstate n, as defined in (j4.2p . The case of the Gibbs ensemble therefore yields the asymptotic relation 

, . En PnTn _ / e~ " " / ^ dpd<j) 

^ ' En^n ~ !e-"^/Tdpd^dX- 

We compare in Remark 18.21 our model of stochastic data with a more standard model having given prob- 
abilities to be in mixed states, which is not an equilibrium solution of the Ehrenfest dynamics. To sample 
from the Gibbs equilibrium density is standard in classical Hamiltonian statistical mechanics but it seems 
non standard for Ehrenfest quantum dynamics. 



8.0.1. The Constrained Stochastic Data. As in models of heat baths [551 HI] and [55] we assume that the 
data of the light particles (here the electrons) are stochastic, sampled from an equilibrium distribution of 
the Liouville equation. All states in this distribution correspond to pure eigenstates of the full Schrodinger 
operator with energy E. There are many such states and here we use the canonical ensemble where the data 
is in state (p with the Gibbs-Boltzmann distribution proportional to exp{—HE/T)d(jfd(f>^dpdX, i.e. in any 
state (f>, for \\<j>\\ = 1, with probability weight 

Let us now determine precise properties of this distribution generated by the Hamiltonian He- To reduce 
the complication of the constraint cj) ■ (f) — 1, we change variables = 0/(0 • 0)^^^ and write the Hamiltonian 
equilibrium density as 

(8.2) exp ( - (p . p/2 + Ao + /t) d^^'d^'dp dX. 

^ • ^ 
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Diagonalize the electron operator V{X'^) by the normahzed eigenvectors and eigenvalues {pj,Xj} 

^■V{X')^ = Xo+J2{Xj-Xo)\1j\^ 

— : Aj 

where 

(8.3) (T/-Ao)pj =XjPj, 

Ao -0, 

with real and imaginary parts 7j =: 7J + 17*. The orthogonal transformation 'l2j"fjPj shows that the 
probability density (|8.2p is given by 



(8.4) V 



using that the determinant of the matrix of eigenvectors is one 



If we neglect the constraint and set J2i>o l7il^ ~ ^ joint distribution density T) is approximated by 



( 








(n,>„/K.e-^.I^.IVTd^;d7j) 





where {7J, 7], J > 0} are independent and each 7J and 7* is normally distributed with mean zero and variance 
T/\j. We see in Lemma |8. II that this approximation is accurate, provided the spectral gap conditions 



(8-5) mmx^GDinax x&d z],>o ' ^ 



Y ^ sX + {\ ~ 3)Xc 
s G [0, 1] 



^J>0 Xj{Y) 



= : K < 1 



hold, where D is the set of attained nuclei positions for the Ehrenfest dynamics. The first condition implies 
that the temperature is small compared to the gap. The second condition means that the electron eigenvalues 
are almost parallel to the ground state eigenvalue, as a function of X. We have dx Aj — pj ■ dx {V—\o)pj which 
typically decays fast with growing eigenvalues, since dx{V — Aq) is smooth and pj become highly oscillatory 
as j — >■ 00. 

Lemma 8.1. Assume the electron eigenvalues have a spectral gap around the ground state eigenvalue Aq, in 
the sense that (j8.5p holds. Then the marginal probability mass 



satisfies 

(8.6) |log4^|<«. 

Proof. We first note that H^illL^ is bounded so that each component |7j| is also bounded. Each integral factor 
has the derivative 

9^/i,„P<c.^-^"'''"''/^'^^^"^"'''^''^rf7;;d7; 



32 



ANDERS SZEPESSY 



for n ^ and the derivative equal to zero for n = 0, so that 
dxr{X) 



' \^) Z^n>0 



'>0 A„(X) /|^„|2<c<="'"'^"'''''''^'^°'^'''''^->^'^->^ 

The integral in the denominator has the estimate 

+0(A„r-ie-^"^/(2^))) 

using eXn/T < |7nP < CXn/T and r/A„ ^ 1 in this domain; the integral over the remaining domain satisfies 

Jl7„P<eA„/T "7ra"7n !i J |2 « «7n "7n 

< (1 + /|^„i2<,A„/(T(i+0) e"l^"l'c^7; ^7^. 

Choose e = 4rA-i log(A„/r) to obtain 

(8.7) / e-^"l-"lV(^i:..o 17.1^)^^. ^7; = f E l7.f (- + 0{T-\-' \og{KT-') 

We have similarly 

/ggN Ji7.P<c T^^^oh.pe - c!7„a7„ 

= £ E,,. J7,f (f + 0(TA-i log(A„r-i))' 



which implies that 



(8.9) |,„s.'W.-./'^- 



/ — (sXc + (1 - s)X) • (X - Xc) ds| < K. 

"'0 



□ 



Remark 8.2. Entropy and the Standard Canonical Density Distribution. Let qj denote the probability of 
electron state j in the data <j). In the setting of the canonical distribution model there holds 

(8.10) q^^er^./T/J2,-h/T^ 

j 

which follows from maximizing the von Neumann entropy defined by — qj log qj , with the probability and 
energy constraints J2j Qj — 1 ^'^.d = constant, see [53]. The stochastic model for the variable |7jp, 

measuring in (|8.4p and (|8.3p the probability to be in electron state j, is different from the model (18.101) . 
since the Gibbs distribution ^^J^h-^itJ^^^^x ~ En PrJnl X^n (|8.3p includes both the density /5„ and the 
weight r„ := |7„p to be in electron state n. 
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8.0.2. The stochastic molecular dynamics models. In this subsection we consider the special case when the 
observable does not depend on the time variable (and the velocity) . Then it is enough to determine an integral 
with respect to the invariant measure; there are several alternatives, cf. 0, and we focus on ensemble averages 
computed by stochastic Langevin and Smoluchowski dynamics. The observable in the Ehrenfest dynamics is 
then 

(8-11) r tY\ -^o(^)/r+iog^S^ 

_ Jr3n g(X)e ■•■(^<:> dX 



_Ao(x)/r+iog-;^ 



where by Lemma |8. II 



|log^|<., 
r{Xc) 



which implies 

J^s.9{X)e-'o(^yMX)dX _ J^,^g{X)e-^oi^y^dX 
^ ' J^,, e-Mx)/T^^x)dX J^s.e-Mx)/Tdx '^^^^'^ 

Let Wt denote the standard Brownian process (at time t) in with independent components and let 
K be any positive parameter. The stochastic Langevin dynamics 

dXt = ptdt 



dpt =-dxXo {Xt)dt - Kptdt + ^/WKdWt 
and the Smoluchowski dynamics 

dXs = -dxX(){Xs)ds + V^dWs 
has the unique invariant probability density 

g^(p.p/2+Ao(X))/T^p^j^ 

respectively 

^~Xo{X)/T^X 

/r3« e"^o(x)/Tdx' 

cf. [9] . We see that the invariant measure for the Smoluchowski dynamics and the marginal invariant measure 
of the Langevin dynamics are equal to the first term in the right hand side of ()8.12p and we have proved 

Theorem 8.3. Both Langevin and Smoluchowski stochastic molecular dynamics approximate Gibbs Ehrenfest 
observables (|8.1ip with error bounded by 0{k), provided the spectral gap condition (|8.5p holds. 

If we assume that 

(8.13) the matrix p„ • GJ2k ^x^ {G^^Pm) is (uniformly in X) bounded in t"^, 

then the difference of the Schrodinger and Ehrenfest Hamiltonians has the bound 

M-i|3?(0 • G Efc Ax, (G- V)) I = M"i|3?( E„ „ lllmPn ■ G A^, (G-^p™)) I 

<0(A/-i)Ej7nP, 
which is asymptotically negligible compared to the Ehrenfest potential energy 

n 

and we obtain 
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Corollary 8.4. Both Langevin and Smoluchowski molecular dynamics approximate Gibbs Schrddinger ob- 
servables 

J g{X)e-"'^''-P'*^^^dXdpd<j>/ J e-'''^^'P''f'y^dXdpdcl> 
with error bounded by 0{M^^ + k), provided (|8.5p and ()8.13p hold. 
By combining (|8.7p and (jS.Sp , the proof in Lemma 18.11 in fact shows 

dx\ogr{X) ^lY.^+O{l)Y,\dxMK'^0gX-' 

n>0 n>0 

which together with ()8.9p and (|8.1ip imply 

Theorem 8.5. If we instead of the gap condition (j8.5p assume the weaker condition 

TA-i _ _ < 1 

(8.14) EnyoldxMK' ^0{1), 

J2n>a\9x,K\K^'^ogXj =K<1, 

then the results in Theorem \8.3\ and Corollary \8.4\ hold, provided the Born-Oppenheimer potential Xq is replaced 
by 

Ao(X) + |^logA„(X). 



2 

n>0 



Note that the correction can be written as 

T 



2 trace-^ log(l/ - Aq) 

where trace""" is the trace in the orthogonal complement of the electron ground state ■00 satisfying V-ipo = AqV'o • 
The work |49| shows that Langevin dynamics, using the rank one friction and diffusion matrix 

K = K{X) = 2M-^/^dxMX) ■ dxMX); 

approximates time-dependent observables based on the Ehrenfest dynamics with the accuracy o(M^^/^) on 
bounded time intervals if k = 0(M~^/^). 

Theorem l8.3l is relevant for the central problem in statistical mechanics to show that Hamiltonian dynamics 
of heavy particles, coupled to a heat bath of many lighter particles with random initial data, can be approx- 
imately described by Langevin's equation, as motivated by the pioneering work [T7],[33] and continued with 
more precise heat bath models, based on harmonic interactions, in [22l[21] [52] . 

9. Discrete spectrum 
This section verifies that the symmetric bilinear form 

vVrV -\- jv'^ dxdXo 



is continuous and coercive on H^{T^'^^^^^'>~^) also for the Coulomb potential, which implies that the spectrum 
of V is discrete by the theory of compact operators, see [18]. Let r := \x^ — Integrate by parts, for any 
e > 0, to obtain 

- Iq W^\v''T-'dr = - Wr'^dr 
^-ifv^dr^dr 



1/2 



J^^vdrv r^dr+r-^YrZo 



> — ^ v'^r'^dr J^^{drv)^r^dr 



[ 2 ir=0 
2 lr=0 
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and integrate the representation v^{R)R — v^{p)R + 2vdrvRdr to estimate the last term 

> /^/, vHr)Rdr - lli^ ( I'^^d^vfr^R'Hr v^R^r-Hr) "'\p 

which shows that 

-Cw^\-'-'dr > ^eC{drv)Wr-y^''v^^r^dr. 
Similar bounds for the other interaction terms in V implies that the bilinear form is coercive 

/t3(j+jv)-i vVrV + dxdXo 

> /t3(.+«)^i I EU I^-^^I' + m Ell |5x"«P + dxdXo, 

for 7 > 60(Af + J). Analogous estimates show that the bilinear form is also continuous, i.e. there is a 
constant C such that 

/ vVrW + jvw dxdXo < C||w|l^i(x3(./+N)-i)|lit;|j/^i(x3(-^+«)-i)- 

The combination of coercivity and continuity in H^{T'^^'^^^^^^) implies, by the theory of compact operators, 
that the spectrum of V consists of eigenvalues with orthogonal eigenvectors in L2^j'3(,/+Af)-i^^ [Tg] . 

Remark 9.1 (The Madelungen equation). An alternative to ()2.7p is to instead include the coupling term 
— 2^ J2j ^ in the eikonal equation, which leads to the so called Madelungen equations [321 ■ Near the 
minima points, where E — Va{X) = 0, the perturbation —-0 • - Aj^j^/; can be negative and then there is 

no real solution dxO to the corresponding eikonal equation. To have a non real velocity dxS is in our case 
not compatible with a classical limit and therefore we avoid the Madelungen formulation. 
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